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Efficient Simultaneous Calibration of a 
Magnetometer and an Accelerometer 

Conrado S. Miranda, Janito V. Ferreira 


Abstract —This paper describes a calibration algorithm to 
simultaneously calibrate a magnetometer and an accelerometer 
without any information besides the sensors readings. Using a 
linear sensor model and maximum likelihood cost, the algorithm 
is able to estimate both sensors’ biases, gains, and covariances, 
besides sensor orientations and Earth’s fields. Results show errors 
of less than 0.1 standard deviations in simulation, and high- 
quality estimates with real sensors even when the algorithm’s 
assumptions are violated. 

Index Terms —Maximum likelihood; Parameter estimation; 
Sensor calibration. 

I. Introduction 

ENSOR reading is an important step in a robot’s control 
loop. However, the values obtained may be inaccurate 
because of incorrect sensor calibration. In the field of aerial 
vehicles, accelerometers and magnetometers are frequently 
used to estimate attitude m. As both these types of sensors 
are calibrated using Earth’s gravitational or magnetic field 
as reference, most research into calibration of these devices 
focuses on magnetometer calibration, which can also be used 
to calibrate accelerometers in stationary settings. This work 
focuses on batch algorithms, that is, algorithms that collect a 
data set and then calibrate the sensors, which can be used to 
provide the initial conditions for online methods that solve the 
calibration problem by filtering o or iterating El, 0. 

A common approach to calibration is to assume that the 
magnetic field’s norm is known or unitary, and then adjust 
the gain and bias so that the readings match this value. 0 
extended the TWOSTEP algorithm, which tries to minimize 
the difference between the field’s norm and the sensed norm, 
to calibrate both the bias and gain. 161 proved that noiseless 
magnetometer readings, including soft- and hard-iron effects, 
occur on the surface of an ellipsoid manifold, and then 
developed an algorithm to calibrate the bias and gain. 

Although both methods can be used to calibrate accelerom¬ 
eters as well as magnetometers, the calibration must be per¬ 
formed independently, so the magnetometer and accelerometer 
are calibrated in different frames, and the rotation between 
them is unknown. Furthermore, neither method is able to 
estimate the direction of the magnetic field, which may be 
unknown and is essential for attitude estimation using magne¬ 
tometers El, 13, 0. 

0 partially solves this issue by calibrating the magne¬ 
tometer while also capturing a reference for the z direction, 
provided by some independent filter. This allows the direction 
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of the Earth’s magnetic field to be estimated, and the cali¬ 
bration to be performed in some previously known reference 
frame, while also estimating the magnetometer’s bias and gain. 
However, a reference orientation, whether reliable or not, is 
not always available, thus limiting the use of algorithm. 

ca tried to solve both problems by simultaneously cal¬ 
ibrating the magnetometer and accelerometer, allowing the 
magnetic field, the attitude, and all the sensors’ parameters to 
be estimated. This eliminates the need for a reference so that 
the calibration can be performed from sensor readings alone, 
that is, no external apparatus or knowledge besides the sensor 
readings is required. However, the problem is replete with low- 
quality local minima, which sometimes significantly degrade 
the calibration performance, and the solution has a very high 
computational cost, making the algorithm impractically slow. 

The novelties in this paper are the improvement over ifTOll 
made possible by the use of better initial estimates of the 
parameters, leading to better minima, and the approximation 
for the rotation estimation, which also helps to avoid poor 
minima and reduce computational time because of its closed- 
form solution, making it suitable for use in practice. The 
initial conditions are given by an adaptation of 0, where 
the z reference is replaced by an estimate computed from the 
accelerometer. To our knowledge, the algorithms described in 
the present work and in ifTOl are the only accelerometer and 
magnetometer calibration algorithms that are able to perform 
the complete calibration from collected data alone. Moreover, 
the algorithm proposed here is also the most statistically 
correct and generic solution based on linear models since all 
the parameters (gains, biases. Earth’s fields, and orientations) 
must be estimated at the same time. 

The algorithm is tested in simulated and real sensors. The 
simulations show that the algorithm correctly calibrates differ¬ 
ent sets of sensor parameters and allows comparison between 
the estimates and the ground truth. Several variations of the 
algorithm are tested, and the best-performing one is also used 
to calibrate the real sensors. The data capture was designed to 
violate the basic assumptions underlying the algorithm so that 
its robustness to adverse real-world conditions could be tested. 
It is shown that even in this imperfect setting the proposed 
algorithm is able to achieve high-quality calibration. 

The paper is organized as follows. Sections HI] and |III| 
describe the sensor model used and the collection and prepro¬ 
cessing of the data, respectively. Based on the sensor model, a 
maximum likelihood cost function is described in Section |iy| 
The initial parameters’ estimates are given in Section |V] 
and the optimization steps are described in Section [Vll The 
results of simulated and real-world experiments are reported 
in Section IVIII and are used to evaluate the proposed method. 
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The concluding remarks are presented in Section I VIII I 

II. Sensor Model 

The linear sensor model with Gaussian noise, where the 
correct value is scaled by a matrix and added to a bias and 
a noise, is frequently used in the literature Q, 0 and is 
also used in this paper. The value Vs G read by the 
sensor s, where s = a and s = m for the accelerometer 
and magnetometer, respectively, is given by; 

Vs= KsV* + bs + es, EsA/'(0, Ss), (1) 

where Kg G and hg G R^ are the gain matrix and bias 

vector associated with the sensor, v* G R^ is the nominal value 
to be read, and Sg G is the noise covariance, which is a 
positive-definite matrix. It is important to note that the linear 
sensor model for the magnetometer is also able to handle soft- 
and hard-iron effects, which can be embedded in the sensor’s 
gain and bias 0. 

Conversely, given the parameters and a value read by the 
sensor, the nominal value can be estimated as: 

vt=K-\vg-bg) + e:, e: 

where it is assumed that the gain matrix Kg is invertible. If 
it is not, then its columns are not linear independent and the 
measured values Vg must lie in a reduced subspace spanned 
by these columns. Therefore, only the case where Kg is 
invertible is considered, and the nominal measurement u* can 
be recovered. 

The nominal value for a sensor is given by the Earth’s 
respective field rotated by the current orientation of the sensor 
frame, so that v* = i?f J, where is the field value for a 
given sensor in the fixed inertial frame I and R G SO{3) is 
the rotation representing the current sensor orientation. Since 
the choice of X is arbitrary, one can choose it such that the 
gravity g only has a negative component in z and the magnetic 
field h has a positive component in x, such that 

5 = [O 0 gz]^, h= [hg; 0 ^, (2a) 

gz <0, K> 0, (2b) 

Sa =9^ Sm = h. (2c) 

It should be noted that the assumption that the magnetic 
field component hx is greater than zero can be satisfied 
almost everywhere except at the magnetic poles, where the 
magnetic field is perpendicular to the Earth’s surface and the 
magnetometer does not provide more orientation information 
than the accelerometer does for the gravitational field. In this 
particular case, one cannot use this combination to perform 
the calibration, but one can do so in all other contexts, as will 
be shown. 

As discussed in cni and highlighted later, the sensor frame 
also has multiple possible definitions, since the gain matrix 
can be written as Kg = RK'g, where i? is a rotation matrix. 
Therefore, some constraints must be fixed to make the solution 
unique. This paper assumes that Ka is triangular, which is 
a sufficient assumption to define a single solution. Moreover, 
Km is allowed to have any value, so that it can incorporate the 
misalignment between the accelerometer and magnetometer. 


Since during the calibration the algorithm does not use 
information about the body frame, and this may not be 
available if the calibration is being performed outside the 
desired body, other algorithms such as mi can be used to 
align the sensor and body frame. 


III. Data Preprocessing 

Assuming the linear model introduced in Sec. |II] and that 
the sensors can be held relatively still when readings are being 
taken so that a set of sensor measurements have the same 
nominal value, the measured values can be written as: 

Vg[i,j]-^N'{g.g[i],^g), g,g\i]= KgRiV^ + bg, (3a) 

jG{l,2,...,A,[*]}, iG{l,2,...,iV}, (3b) 


where N is the number of measurement sets, As[f] is the 
number of data points collected for the i-th set, g,g[i] is the 
mean value shared by the measurements in the i-th set, and 
Ri is the orientation at which the i-th set was collected. 

Since the measurements in each set are assumed to have 
sampled the same information, the mean values gg[i], which 
are different for each set, and the covariance matrix Eg, which 
is shared among all sets, can be estimated. 

The minimum variance unbiased estimator for g,s[i] is the 
sample mean ifT^ . given by: 


fis\i] 


1 

Ag[i] 


A,[i] 

i=i 


(4) 


and has a distribution described by: 

(5) 

Since the sample mean jlg\i\ characterizes the values for the 
2 -th set, it can be used instead of them as a single measure¬ 
ment with appropriate weight, associated with the number 
of measurements As[ 2 ], during optimization to decrease the 
computing time. 

The minimum variance unbiased estimator for the noise 
covariance is given by; 

^ Eti - A. W) j] - 9s[i]f 

— / M \ 5 

where the subtraction of N from the denominator is a gener¬ 
alization of the Bessel’s correction m for N measurement 
sets, instead of 1 as in the original correction. 


IV. Cost Eunction 

Given that it is assumed that only the measured values are 
known, the full parameter set is given by all the sensors’ and 
fields’ parameters and all orientations. This set can be written 
as; 

0 Kd , ba , gz; 7 Km 7 bm 7 hx 7 hz 7 Rl , R2 7 7 }• 

Note that this is the most generic setting possible for this 
problem with the assumption of affine measurements, since 
all parameters are considered unknown and any previously 
known parameter can be fixed to the known value. 
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Using the model defined in Eq. Q, the negative log- 
likelihood can be used to create a cost function m, which is 
given by: 


By separating the known terms in Eq. ( l9ab from the un¬ 
known terms and considering all N sample sets, the approxi¬ 
mation can be written as: 


0 ^( 0 , 5 ^, 5771 ) — 

N A«[i] 

s€{a,m} 2=1 i=l 

Vs [i,j] = KsRiVs +bs-Vs [ij] ■ 


(7a) 

(7b) 


In this format, the calibration problem becomes a traditional 
minimization problem over all the parameters in 0 , where the 
restrictions over g^, h^, and Ri must be satisfied. 

A simplified cost function can be defined on a smaller 
parameter set 0 ' = 0 \ {Eq, E^}, where the covariances are 
replaced by their estimates computed using Eq. (| 6 ]l. In this 
case, the maximum likelihood cost can be simplified to 


N 


J'(0';/ia,/im,So,E,„) = ^ ^ As[i]/is[z]'^E 3 

2=1 s€{a,m} 

( 8 a) 


p.s[i] = KsRiV^ + bs- fis[i\, 


( 8 b) 


J 77 = 




Jn 


77 0 


Ji = [AsW^ ® As[*]^ 1] 

vecA 


77 = 


( 10 ) 


where < 8 ) denotes the Kronecker product and vec denotes the 
vectorization operator. Assuming that A is symmetric, the 
nontrivial solution 77 that minimizes the approximation error 
is given by the right eigenvector corresponding to the smallest 
singular value of J 191. 

Since any arj is also a solution to the approximation in 
Eq. (O, one can find the correct a by manipulating the 
elements in Eq. (|9ll, which provides the correct scale. Hence 
the correct value is given by: 


in which the real measurements Vs[i:j] can be replaced by 
their sample means As[*] with weight As[ 7 ], which greatly 
reduces the computational cost. 

Since both cost functions are subject to many local minima, 
which can lead to poor estimates as shown in ITOl . a good 
initial estimate must be provided. Then, instead of optimizing 
all the parameters at once, it will be shown that optimizing 
subsets one at a time greatly simplifies the problem. The next 
two sections deal with these problems. 


a = 

Erom the definitions in Eq. (|9]l and the estimates 77 and a, 
the gain and bias must satisfy 

K-^K-^=aA ( 11 ) 

bs = 

S 2 ’ 


V. Initial Estimate 

This section is an adaptation of the work presented in 0, 
where a magnetometer was calibrated using a known correct 
reference measure of the z direction. Since it is assumed that 
no such z is available, this reference is replaced here by 
the estimated accelerometer’s field value, which only has a 
negative component in the z direction. 

Eurthermore, the constraint that the magnetometer cannot 
have mirrored axis relative to the inertial axis is relaxed, 
which was not possible in the original reference. The reader 
is referred to the original paper for a more detailed derivation 
of the equations and implementation details. 

Assuming that the field is unitary and that the noise 
eg can be neglected, the following approximation can be 
performed: 

0 Ri As W ^^As W + As [*] + c (9a) 

A = K-^K-^ (9b) 

b=-2bjK-^K-^ (9c) 

c = b'^K-^K;^bs-l, (9d) 

which follows directly from taking the norm of both sides of 
Eq. ([TJ and replacing the single measurements by their mean 
given by Eq. (|4|, which has smaller variance and thus is a 
better estimate. 


which gives the bias estimate for each sensor directly. 

The gain Kg can not be uniquely determined, as any KgR 
with RR^ — I 3 is also a solution to Eq. (fTTT i. If the upper 
triangular solution obtained by the Cholesky decomposition is 
denoted by Kg, the estimates can be written as: 

Km = KmR, Ka = Ka ( 12 ) 

where it is assumed that only Km is subject to rotation as the 
rotation between the two gains is fixed and there would be 
multiple solutions if both had full unknown rotations. 

So far, the solution described in 0 , which provides esti¬ 
mates for the sensor bias and gain, has been used to compute 
the initial estimates. This solution could be applied here 
since the accelerometer and magnetometer could be consid¬ 
ered decoupled problems. However, in Eq. (fT2l i. Ka can be 
estimated directly while Km depends on the unknown rotation 
R. Since only the magnetometer was calibrated in 0 , this was 
not an issue as the ground-truth z direction was available. 
Moreover, the magnetic field component hz could also be 
estimated directly. To solve this problem, the ground truth z 
was replaced by the estimate of the nominal value v* from 
the accelerometer, since this is already calibrated. 

Therefore, estimates for the rotation R between the ac¬ 
celerometer and the magnetometer and the component of 
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the magnetic field can be found simultaneously by solving the 
following optimization problem: 

1 ^ 

min {Aa[i\ + Amii]) \\hz + Za\i]'^ZmiilWl (13a) 

’ * i —1 

s.t. Zs[i] = k~'^ (13b) 

R e 50(3) (13c) 

where and Am are used to weight the number of samples 
in each data set. As stated earlier, this equation is similar to the 
one in 0, but the ground truth z is replaced by the estimates 

Zal'i]- 

Test showed that, unlike in 0, this problem is subject to 
a few low-quality local minima because of the replacement 
of precise knowledge of the z direction by the accelerometer 
estimate and that the initial condition R — I 3 and hz — 0 
sometimes leads to poor estimates. Hence, an optimization 
with random restarts was performed and the initial conditions 
were given by a random rotation as described in na and 
hz ~ W([—1,1]). Tests showed that at most 100 iterations 
were needed to find good initial estimates. 

Since it was assumed that the fields were unitary, it follows 
that 

k = 

As discussed in na, having the fields’ and gains’ scales as 
variables leads to multiple solutions, since the gain can be 
multiplied by a constant while the field is divided by the same 
value. Therefore, in this paper the value Pz = — 1 and hx = ^ 
are adopted, and the magnetometer estimates must be adjusted 
as: 

km ~ kmhx-, h^ — hz/hxi 

while the accelerometer does not have to be adjusted because 
it already has a unitary field. 

Finally, the rotation Ri associated with each sample set 
is given by the eigenvector corresponding to the minimal 
eigenvalue of Bi M, where Bi is defined as: 

Ri — 2 (14a) 


= f{K, 

■ + 

4 - 

S^s) 


(14b) 

= K7^ 

ihs [i] 

-bs), 

Ws = 


(14c) 


'0 

-2/1 

-2/2 

-2/3 



fix,y) = 

yi 

0 

-2:3 

r\ 

X 2 


(14d) 

y2 

X 3 

u 

-Xl 




ys 

-X2 

Xl 

0 




and the weights Wa and Wm introduce uncertainty weighting. 

All initial estimates therefore have closed-form efficient 
solutions except for R and /i^,. However, Eq. (fOT l has few 
decision variables and can be efficiently optimized by random 
restarts and gradient descent. 

VI. Optimization Steps 

As the cost functions in Eqs. (I7]i and (O have many joint 
parameters, the parameter set is split into four disjoint subsets 
that can be solved efficiently one at a time. The optimization 
starts with the estimates provided in Sec. |V] and then iterates 


each subset optimization while keeping the other parameters 
fixed until convergence is achieved. Convergence is defined by 
Jk-i — Jk < 1, where Jk is the cost after the fc-th iteration 
and 7 is the stop parameter. Since the true values for any 
parameter in 0 are not available, the hat notation is dropped 
where convenient to simplify the expressions. 

This section is an extension of the optimization presented 
in Col. The main differences are, firstly, the approximation 
for the rotation estimation described in Sec. IVI-AI This uses 
Eq. (fl4l) to compute an efficient approximation to the real ori¬ 
entation during the optimization, which is not present in ifTOI . 
leading to a significant slow-down of the algorithm and worse 
results. Secondly, the format of the gain matrix Ka computed 
in Sec. IVTCl is assumed to be triangular here because its initial 
condition is given by a Cholesky decomposition while in ifTOll 
it is assumed to be symmetric. The main difference between 
the method proposed here and the one described in 0, apart 
from the fact that two sensors are used and no ground-truth is 
available in the former, is that the optimization is performed in 
steps with efficient solutions, while in 0 the initial estimate is 
computed and then all parameters are optimized at the same 
time, preventing the inherent structure of the problem from 
being exploited. 

A. Rotation Estimation 

Erom Eq. ((Sj), it is clear that the cost for the z-th rotation is 
given by: 

4, = ^ (15) 

/is[z] = KsRiS^ + bs- As[*]- 

As a closed-form solution could not be found, the problem is 
optimized directly through gradient descent using a quaternion 
to represent the rotation. However, the approximation used to 
compute the initial estimates, which is given by the eigenvector 
associated with the minimum eigenvalue of Bi from Eq. (fT4l) . 
can be used to speed-up the initial steps of the optimization, 
when the estimates are still far from the true values. The 
performance of this approximation approximation is evaluated 
in the experimental section. 


B. Bias and Field Estimation 


Since the sensor readings are assumed linear in the biases 
and fields and have Gaussian noise, the problem can be posed 
as a generalized least squares (GLS) problem El given by 
finding (3 to approximate 

2 / = A/3 -f e, E[e|A] = 0, Var[e|A] = H, 


whose solution is given by: 

/3= (A'^H-^A)”^ 

Writing jla[i\ in the GLS format using Eqs. 0 and 0 gives 
KaRlZ I 3 


Aa[l] 


Mn 



H = diag 


R-aRN Z 13^ 


Aa]!]’’"’ Aa[A] 
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where z is the unitary vector in the z direction. As Sec. [Vl 
established that = — 1 , the gain matrix is adjusted to 
ensure that this equality is still valid. The estimate is therefore 
adjusted so that K'^ = —K^gz, where g^ is the estimated 
gravity. 

The magnetometer estimate fim [i] can be written in a similar 
fashion, allowing h^, hz and bm to be estimated. As = 1, it 
follows that K'^ = Krnhx and h'^ = hz/hx, which maintains 
the scale. 


C. Gain Estimation 


The gains can also be written in the GLS format using 
Eqs. Q and (|5]l. Assuming that Ka is upper triangular, which 
is valid for the initial estimate since it comes from a Cholesky 
decomposition, one has that: 




'Gi 


X.n' 

Ka,12 




Ka, 13 

_fLa[N] - ba_ 


_Gn_ 


Ka,22 

Ka,23 

Ka,33 


G, = 


9Ri,i 

0 

0 


9Ri.2 

0 

0 


9Ri,3 

0 

0 


9R,] = {R 9 )], ^ = diag 


9Ri,3 

0 


0 

0 

9Ri,3 


9Ri,2 

0 

Sa 

A^TV] 




where {v)i represents the *-th component of the vector v. 

Again, a similar formulation can be used for the magne¬ 
tometer, but the gain Km is considered full, as it is the product 
of a rotation with an upper triangular matrix. 


D. Covariance Estimation 


The maximum likelihood covariance estimator is similar 
that given in Eq. ([ 6 ll, but the reference is replaced by 
an estimate of g,s\i] using parameters in 0, and Bessel’s 
correction is not used, since it lowers the likelihood of the data 
in order to provide an unbiased estimate. The new estimate is 
then given by: 


S. = ■ 

Ss[iJ] 


= Ss[iJ]- {KsR^s^ + bs) . 


(16) 


The performance with test data, which is not used for 
training, can be evaluated to test the generalization of the 
algorithm iflTl . 

In the real experiment, a magnetometer and an accelerom¬ 
eter had their readings measured while being held by hand 
inside a building, which slightly violates the hypothesis that 
all samples measure the same value. It will be shown that the 
algorithm is able to fit the sampled data correctly. 

It is important to note that only the algorithm derived 
in this paper is used for calibration since the authors do 
not know of any other work, apart from their own previous 
contribution Eol upon which this paper is based, that describes 
an algorithm that can perform this kind of calibration without 
external references. If an algorithm such as the one proposed 
in 13 were used, additional information would have to be 
provided, making comparison of the results impossible since 
the calibration would be performed on different data. Thus, no 
comparison with any previous work is performed. 


A. Simulated calibration 

To evaluate the calibration algorithm, 100 Monte Carlo 
simulations with new sensor parameters on each run were 
performed, allowing evaluation of both common and possible 
exceptional calibration behavior. This analysis is important to 
show that the method described handles different sensor char¬ 
acteristics well. Table U presents a summary of the parameters 
used in the experiment; the symbols in the table are those used 
throughout the paper. 

Eor each Monte Carlo run, a total of A = 15 sets of 
sensor measurements were considered, with the same number 
of samples A[i] ~ ^^({400,401,..., 600}) for both sensors. 
The rotations for each set were randomly created using the 
method described in Ea, and the stop condition was set to 
7 = 10-1 

When varying the stop condition 7 , the same sample sets 
were used for each value tested, so that the performance can 
be viewed as the same algorithm stopping at different times. 
This avoids differences caused by random fluctuations and 
allows better comparison. Eor a varying number of sets, runs 
with N + 1 sample sets merged a new set to the previous N 
already used, simulating collection of an increasing number of 
sets and also preventing random fluctuations from significantly 
affecting the comparison. Alternatively, the training with N 
sets can be viewed as training on a subset of the data set 


VII. Experiments 

Two experiments are performed to validate the proposed 
algorithm: calibration of simulated sensors and calibration of 
real sensors. The simulation is used to evaluate the effects of 
changing the stop condition, the number of sampled intervals, 
and the algorithm hypothesis. The performance is analyzed in 
terms of errors in relation to the ground-truth values, which 
are available in the simulation, and the total time required 
for the simulation. Using a Monte Carlo approach, many 
combinations of sensors can be simulated to evaluate the 
robustness of the proposed algorithm. 


Table I: Experimental parameters 



Description 

Value 

N 

# of measurement sets 

15 

A[il 

# of samples for each set i 

W({400, 401, ...,600}) 

Ri 

Rotation for set i 

See (B] 

7 

Stop condition 

lo-'^ 

9z 

Gravity’s 2 component 

W([-1.5,-0.5]) 


Magnetic field’s x component 

W([0.5,1.5]) 

hz 

Magnetic field’s 2 component 

W([-1.5,1.5]) 

^s,i 

Bias’ component 

w{[-i,il) 

Ks 

Gain mati'ix 

See Eq. {\1\ 

Ss 

Covariance matrix 

See Eq. G8J 


# of Monte Carlo runs 

100 
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Figure 1; Reconstruction error for training data for different stop conditions. Best viewed in color. 
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Figure 2: Reconstruction error for test data for different stop conditions. Best viewed in color. 


provided to the algorithm that trained with N + 1 sets, thus 
maintaining consistency. 

The field components gz, hx and hz are sampled uniformly 
over [—1.5,—0.5], [0.5,1.5] and [—1.5,1.5], respectively. Pa¬ 
rameters for the accelerometer and magnetometer are created 
in a similar fashion. Each component 6^,1 of the biases has a 
random value sampled uniformly in [—1,1]. The gain matrices 
are given by 

K, = l3+ SK,, SK, ~ Z^([-0.1, 0.1]), (17) 

and the covariances are given by 

if ([0.5,2]), ifz=j 
|Eg_ij ^ if ([— 0 . 2 , 0 . 2 ]), otherwise, 

where as = 10 *^, e ~ if([— 2 ,—4]), which guarantees that 
the covariance matrices are positive definite, but allows high 


condition numbers. Instead of using the nominal value K^, 
the magnetometer is rotated by R' and mirrored by M, so 
that the resulting gain matrix is given by = MR'Km, 
showing that the algorithm described is able to handle any 
configuration between the accelerometer and magnetometer, 
unlike Q. The rotation R' is chosen randomly 1(151 . and M = 
diag(5i, ( 52 , ^ 3 ), 5i ^ ff({—l, 1 }), is a random diagonal matrix 
that indicates the mirroring of each component. 

After calibration, the estimated value of 0 was used to 
rebuild estimates /is[z] for /is[*], as in Eq. (O. To compare 
the quality of the reconstructions, the following measurement 
was used; 


5.4 = 


\ 


(FsH - Eg ^ {fis\i] - 

EtiAW 


( 19 ) 
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Figure 3; Reconstruction error for training data for different numbers of sample sets. Best viewed in color. 
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Figure 4: Reconstruction error for test data for different numbers of sample sets. Best viewed in color. 


where Eg is the real sensor covariance. This measurement is 
derived from the Mahalanobis distance lITSl and represents the 
average error in standard deviations. 

The rotation optimization described in Sec. IVI-AI can be 
performed in one of two ways: always using gradient de¬ 
scent in Eq. (fTSl l or using the approximation in Eq. (fTTl i 
until Jfe > Jk-i and then changing to gradient descent in 
subsequent iterations. As the approximation does not minimize 
the cost directly, eventually the change to gradient descent 
will be made. Another variation involves re-estimating the 
covariance using Eq. ( fTbl l or only using the initial estimate 
given by Eq. (| 6 ]l. As some approaches consider the covariance 
matrix to be diagonal [| 6 l, this restriction is also evaluated. 

Hence there are four algorithms to compare: not refitting 
covariance and optimizing the rotation directly using gradient 
descent (NCDR); not refitting covariance and approximating 


the rotation (NCAR); and fitting full covariance and diagonal 
covariance and approximating the rotation (ECAR and DCAR, 
respectively). The results obtained while readjusting the co- 
variance with direct rotation optimization are not reported as 
these were considerably worse. 

Eor each sample set, another set with different samples and 
rotations for the same sensor parameters was created to eval¬ 
uate the generalization of the fitted parameters, that is, if the 
parameters are able to fit the data not present during calibration 
well. These new samples follow the description in Sec. |III] 
and provide new values which are used to compute the 
new rotations. The rotations were first approximated and then 
optimized directly to reduce the cost function, as described in 
Sec. IVTAl 

Eigure [T] shows the reconstruction error for different stop 
conditions. Eor 7 < 10“^, the error does not change signifi- 
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candy, justifying the choice of this value as a reference when 
evaluating different number of samples. Except for NCDR, all 
the algorithms have similar behavior for lower values of 7, 
while adjusting the covariance was beneficial for larger stop 
conditions. However, as larger stop conditions can only be 
justihed if the computational cost is very high, it will be shown 
that this difference is not relevant. NCDR has slightly worse 
magnetometer calibration because of the failure of the gradient 
descent method to escape local minima. 

It is interesting to note the apparent magnetometer over¬ 
fitting for 7 < 10“^, as the reconstruction error for the 
training set increases even though the cost function decreases. 
Nonetheless, Fig. |2] shows that the reconstruction error for the 
new data not used during training generally decreases. This 
suggests that, although the reconstruction worsens, the param¬ 
eters improve the generalization. The signihcantly higher error 
for one case with NCDR indicates that it is not as robust as the 
other methods using approximations. For all the algorithms, 
about 75% of all cases have errors of less than 0.1 standard 
deviations from the correct value for both sensors, showing 
that the proposed method is able to fit the parameters well. 

Figure [2 shows that all the algorithms eventually converge 
to solutions that are close to each other as the number of 
samples increases, with errors on every run lower than 0.1 
standard deviations for TV > 25 for all methods, except for 
one run with NCDR. However, a diminishing return on the 
number of samples is observed for N > 20 as the error 
decreases very little if more data is collected. The NCAR 
algorithm appears to be the most robust, with a lower error 
spread for a smaller number of samples, so that less user 
data is needed to achieve acceptable performance for most 
situations. The NCDR algorithm usually has the largest errors, 
which, as in the previous case, are the result of poor local 
minima. Both algorithms with readjusted covariance have 
similar performance, indicating that the diagonal simplihcation 
is a good one. Figure |4] shows slightly better performance 


overall, indicating that the estimates obtained are able to 
generalize and ht new unseen samples well. 

While the NCAR, FCAR and DCAR algorithms had similar 
performances for reasonable stop conditions and numbers 
of samples, the NCDR algorithm not only had the worst 
performance for all the conditions simulated but also required 
more computational time, as shown in Fig. |5a| It should be 
noted that, although DCAR has a simpler covariance to ht 
than FCAR, it has a higher computational cost because of 
implementation details. Hence, the calibration time should be 
about the same for these algorithms in a hnal implementation. 
As stated earlier, using higher values of 7 and covariance re¬ 
estimation would be justihed if the computational cost was 
high. However, the cost for NCAR is signihcantly lower, with 
7 = 10“^ having the same cost as 7 = 10 for FCAR. 
Therefore re-estimating the covariance using the rotation ap¬ 
proximation does not improve performance and instead slows 
the processing considerably. 

As shown in Figs. [T] and |2] 7 values of less than 10“"* 
do not improve the performance signihcantly but cause the 
computing time for the approximation methods to increase. 
This occurs because rotation optimization is still taking place, 
which, although having little effect on the calibration, is time 
consuming. This also explains why the time taken by the 
NCDR algorithm stabilizes. 

The time for NCAR is almost constant for all the numbers 
of samples tested, as shown in Fig. |5b] This occurs despite the 
fact that the cost function in Eq. ® is a sum that increases 
as the number of sample sets increases and the stop condition 
compares two iterations of the algorithm. As more terms are 
used in the sum, the difference Jk-i — Jk would also be 
expected to increase, so that it takes longer to satisfy the 
stop condition. Moreover, the optimization steps also increase 
in complexity, requiring more time for each step in Sec. [Vl] 
Nonetheless, the presence of more samples increases the speed 
at which the estimates converge, requiring less iterations to 
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(a) Accelerometer (b) Magnetometer 


Figure 6: Calibrated real sensors. The light lines show the values read and the thick, darker straight lines show the estimated 
reading for each set. Best viewed in color. 


achieve smaller values of J. 

B. Real calibration 

Once it has been shown that the proposed algorithm works 
correctly with different simulated sensors, calibration of real 
sensors must be tested. The main difference between the real 
and simulation experiments is that, although the simulation 
allowed a wide variety of sensor parameters to be tested, it 
also simulated a perfectly still system during the capture of 
each dataset, which was an assumption made in Sec. |III] Since 
the NCAR algorithm showed the best performance across all 
simulations, it was chosen to calibrate the real system. 

In real-world applications the system may not be perfectly 
still, and the robustness of the algorithm in applications where 
this is the case must be tested. Another assumption was that 
the magnetic field is constant throughout the measurements, 
since it is assumed that it can be written as in Eq. (|2]i. This also 
may not be valid in the real world because of the influence of 
power cables and any metal structures nearby. 

To test the algorithm’s robustness under these conditions, 
a RoboVero board was held by hand to introduce vibration 
and measurements were taken inside a large room so that 
magnetic disturbances were present but not overwhelming. 
Both the accelerometer and the magnetometer were configured 
to have the same sample rate so that both sensors took the 
same number of measurements for each set. The set length 
was determined by using time slices in which the norms of 
the values read were kept approximately constant. This also 
introduces errors as the underlying values can vary without 
changing the norm considerably. 

The results are shown in Fig. |6] where the darker lines 
are the estimated mean sensor readings for each measurement 
set. These are computed from the estimated parameters using 
Eq. ®, with length equal to the number of measurements 
As[i] in each set. The light lines are the real sampled values. 


Fig.|6a]shows that the estimated values for the accelerometer 
were very close to the real ones despite the sensors not being 
held still during each set. This confirms the robustness of the 
algorithm when the first assumption is violated and supports 
its use when the sensors cannot be kept stationary. 

Fig. |6b] shows that the magnetometer is also calibrated 
correctly but that its measurements vary a lot more than the 
accelerometer measurements. As a result, some measurement 
sets have slight offsets, as is the case of the sets around 20000 
and 26000 samples. These errors are associated with distortion 
of the magnetic held by the building’s metallic structure. 
However, the small size of these offsets show that this does not 
affect significantly the calibration, confirming the algorithm’s 
robustness in the presence of this kind of distortion. 

It will be recalled that the calibration was performed using 
only the measured values shown in Fig. |6] without any external 
references or sensors. Therefore, the proposed algorithm’s 
robustness to violation of the assumptions about the only mea¬ 
surements available, together with the high-quality calibration 
that it provides, should make it a good candidate for any 
calibration of these kinds of sensors. 

VIII. Conclusion 

This paper has described an algorithm to calibrate the 
parameters of an accelerometer and a magnetometer using only 
sensor measurements without any external information. The 
algorithm is designed to estimate the gain, bias, and covariance 
of each sensor, as well as the orientation of each measurement 
and the direction of the Earth’s magnetic held. Hence, the 
calibration can be performed almost anywhere by anyone as 
only the sensor readings are required. This is the most generic 
setting for this kind of problem, since all the parameters are 
calibrated and any external information can be considered a 
constraint on the parameters and does not make the calibration 
harder. 
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A comparison was made of the base method and its variants, 
which use an approximation of the cost function being mini¬ 
mized. Simulated results show that the simplest variant is also 
the fastest and most robust, with the smallest worst-case errors. 
All the algorithms are able to achieve an eiTor of less than 0.1 
standard deviations for the data used to calibrate the devices, 
indicating that the algorithms were correctly trained, and also 
for new data not used in the calibration, indicating that the 
parameters obtained are suitable for use when taking sensor 
readings and that the learning algorithm was able to generalize. 
Test with real sensors and adverse conditions that violate 
the assumptions underlying the derivation of the algorithms 
showed that high-quality calibration can be achieved in non- 
controlled settings, confirming the robustness of the algorithm. 

Future studies should investigate use of the estimated sensor 
rotations to compute the rotation between body and sensor 
frame, as described in ca, since the sensor frame used was 
specific to this study and may be different from the desired 
frame. Another possibility for further study is to integrate 
external references, such as vision systems m, so that other 
sources of information can be used during calibration to reduce 
the errors even further. 
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